rm(list = ls())
require(mfx)
require(zTree)
require(stargazer)
require(stats4)
require(lmtest)
require(bbmle)
require("multiwayvcov")
require(msm)
require(nlWaldTest)
require(sandwich)
require(miceadds)


zTT1= zTreeTables("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/Ambueletalexperiments/Data/October20ChangingPriors.xls")
#zTT1= zTreeTables("C:/Users/nneligh/OneDrive - University of Tennessee/Dean Stuff/Ambueletalexperiments/Data/October20ChangingPriors.xls")
subjects1=zTT1$subjects

zTT2= zTreeTables("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/Ambueletalexperiments/Data/October23ChangingPriors.xls")
#zTT2= zTreeTables("C:/Users/nneligh/OneDrive - University of Tennessee/Dean Stuff/Ambueletalexperiments/Data/October23ChangingPriors.xls")
subjects2=zTT2$subjects
subjects2$Subject=subjects2$Subject+100

zTT3= zTreeTables("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/Ambueletalexperiments/Data/October24ChangingPriors.xls")
#zTT3= zTreeTables("C:/Users/nneligh/OneDrive - University of Tennessee/Dean Stuff/Ambueletalexperiments/Data/October23ChangingPriors.xls")
subjects3=zTT3$subjects
subjects3$Subject=subjects3$Subject+200

subjects=rbind(subjects1,subjects2,subjects3)

priors=c(0.5,0.85)

#Clean data
subjects$responder=(subjects$"response1[50]">0)
subjects=subset(subjects,responder>0)

uidlist=unique(subjects$Subject)

#find N
#**
length(uidlist)
#**


#create a frame with observations
obsframe=data.frame(uid=vector('numeric'),block=vector('numeric'), priorA=vector('numeric'), number=vector('numeric'), state=vector('numeric'), response=vector('numeric'))
for(i in 1:length(uidlist)){
playframe=subjects[i,]
uidplace=uidlist[i]

#block 1
for (j in 1:50){
blockplace=1
priorplace=priors[playframe[1,682]]
numberplace=j
stateplace=playframe[1,(735+j)]
responseplace=playframe[1,(71+j)]
addframe=data.frame(uid=uidplace,block=blockplace, priorA=priorplace, number=numberplace, state=stateplace, response=responseplace)
obsframe=rbind(obsframe,addframe)
}

#block 2
for (j in 1:50){i
blockplace=2
priorplace=priors[playframe[1,683]]
numberplace=j
stateplace=playframe[1,(785+j)]
responseplace=playframe[1,(221+j)]
addframe=data.frame(uid=uidplace,block=blockplace, priorA=priorplace, number=numberplace, state=stateplace, response=responseplace)
obsframe=rbind(obsframe,addframe)

}

}

expframe=obsframe

#definitions
expframe$priorB=1-expframe$priorA
expframe$Aresponse=(expframe$response==1)
expframe$Bresponse=(expframe$response==2)
expframe$stateA=(expframe$state==1)
expframe$stateB=(expframe$state==2)
expframe$correct=(expframe$response==expframe$state)

write.table(expframe, "C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/Revision/Reports/AmbhuelChangingPriors.csv", sep = ",", col.names = T, append = F)


#Begin constructing first table
pAgiven2=vector('numeric')
restrictionpAgiven1=vector('numeric')
pAgiven1=vector('numeric')
prob=vector('numeric')

#Base Regression
modelCP=lm(Aresponse~factor(stateA+priorA)+0,data=expframe)

#Adjust COV
vcovmodCP=vcovHC(modelCP, type="HC3")
vcovmodCPcluster=cluster.vcov(modelCP, expframe$uid)

modelcpsterr=vector("numeric")
#Model CP sterrors
for(i in 1:4){
modelcpsterr[i]=sqrt(vcovmodCPcluster[i,i])
}


logitCP=glm.cluster(data=expframe, Aresponse~factor(stateA+priorA)+0, cluster="uid", weights=NULL, subset=NULL,  family="binomial" )

#check for violations
NIASindiframe=data.frame(uid=vector('numeric'),pval0.5=vector('numeric'),pval0.6=vector('numeric'),pval0.75=vector('numeric'),pval0.85=vector('numeric'))
for(i in 1:length(uidlist)){
playframe=subset(expframe,uid==uidlist[i])
frame50=subset(playframe,priorA==0.5)
output50=ols(Aresponse~stateA,data=frame50,robust=TRUE, cluster=NULL)
model50=output50[[1]]
vcov50=output50[[2]]
pAgiven2[1]=model50[1]
restrictionpAgiven1[1]=((2*priors[1]-1)/priors[1])+pAgiven2[1]*((1-priors[1])/priors[1])
pAgiven1[1]=model50[1]+model50[2]
sterr=as.numeric(deltamethod(~x1+x2-((2*0.5-1)/0.5)+x1*((1-0.5)/0.5),model50[1:2],vcov50))
tstat=(pAgiven1[1]-restrictionpAgiven1[1])/sterr
prob[1]=pnorm(tstat)


frame85=subset(playframe,priorA==0.85)
output85=ols(Aresponse~stateA,data=frame85,robust=TRUE)
model85=output85[[1]]
vcov85=output85[[2]]
pAgiven2[2]=model85[1]
restrictionpAgiven1[2]=((2*priors[2]-1)/priors[2])+pAgiven2[2]*((1-priors[2])/priors[2])
pAgiven1[2]=model85[1]+model85[2]
sterr=as.numeric(deltamethod(~x1+x2-((2*0.85-1)/0.85)+x1*((1-0.85)/0.85),model85[1:2],vcov85))
tstat=(pAgiven1[2]-restrictionpAgiven1[2])/sterr
prob[2]=pnorm(tstat)

addframe=data.frame(uid=uidlist[i],pval0.5=prob[1],pval0.85=prob[2])
NIASindiframe=rbind(NIASindiframe,addframe)
}

#need to clean an NaN from first part
cleanpval0.5= NIASindiframe$pval0.5[!is.na(NIASindiframe$pval0.5)]

unexpectedchangers0.5=100*sum(cleanpval0.5<0.05)/length(NIASindiframe$uid)
unexpectedchangers0.85=100*sum(NIASindiframe$pval0.85<0.05)/length(NIASindiframe$uid)

violateframe=data.frame(prior=priors,persigvio=c(unexpectedchangers0.5,unexpectedchangers0.85))
violateframe=t(violateframe)

#***********************************
stargazer(violateframe,summary=FALSE,colnames=FALSE)
#************************************

#Find thresholds and look for dropping

priorchangeframe=data.frame(uid=vector('numeric'),pA0.5=vector('numeric'),pA0.5lower=vector('numeric'),pA0.5upper=vector('numeric'),pA0.85=vector('numeric'),ChangePost0.5_0.85=vector('numeric'),Chisq0.5_0.85=vector('numeric'),Prob0.5_0.85=vector('numeric'))

for(i in 1:length(uidlist)){
playframe=subset(expframe,uid==uidlist[i])
model=lm(Aresponse~factor(stateA+priorA)+0,data=playframe)
na1=model$coef[1]
na2=model$coef[2]
a1=model$coef[3]
a2=model$coef[4]


if(sum(playframe$Aresponse*(playframe$priorA==0.85))==0){
postAhigh0.85=0
}else{
postAhigh0.85=(0.85*a2)/(0.85*a2+0.15*na2)
}

if(sum(playframe$Aresponse*(playframe$priorA==0.5))==0){

postAhigh0.5=0
postAhigh0.5upper=0
postAhigh0.5lower=0

diff0.5_0.85=postAhigh0.5-postAhigh0.85
chisq0.5_0.85=1e7*(postAhigh0.5-postAhigh0.85>0)
diff0.5_0.85p=7.888609e-31*(postAhigh0.5-postAhigh0.85>0)+(1-7.888609e-31)*(postAhigh0.5-postAhigh0.85==0)


}else{
postAhigh0.5=(0.5*a1)/(0.5*a1+0.5*na1)

vcovmod=vcovHC(model, type="HC3")

threshsterr=as.numeric(deltamethod(~(0.5*x1)/(0.5*x1+0.5*x2),coef(model)[c(3,1)],vcovmod[c(3,1),c(3,1)]))
postAhigh0.5upper=postAhigh0.5+1.96*threshsterr
postAhigh0.5lower=postAhigh0.5-1.96*threshsterr

if(sum(playframe$Aresponse)==length(playframe$Aresponse)|sum(playframe$Aresponse)==0){

diff0.5_0.85=0
chisq0.5_0.85=0
diff0.5_0.85p=1

}else{

diff0.5_0.85=postAhigh0.5-postAhigh0.85

if(diff0.5_0.85==0){
chisq0.5_0.85=0
diff0.5_0.85p=1
}else{
test0.5_0.85=nlWaldtest(model,"(0.5*b[3])/(0.5*b[3]+0.5*b[1])-(0.85*b[4])/(0.85*b[4]+0.15*b[2])",Vcov=vcovmod)
chisq0.5_0.85=as.numeric(test0.5_0.85[1])
diff0.5_0.85p=as.numeric(test0.5_0.85[4])
}
}

}
addframe=data.frame(uid=uidlist[i],pA0.5=postAhigh0.5,pA0.5lower=postAhigh0.5lower,pA0.5upper=postAhigh0.5upper,pA0.85=postAhigh0.85,ChangePost0.5_0.85=diff0.5_0.85,Chisq0.5_0.85=chisq0.5_0.85,Prob0.5_0.85=diff0.5_0.85p)


priorchangeframe=rbind(priorchangeframe,addframe)
}

percent=vector('numeric')
percent[1]=100*sum((priorchangeframe$pA0.5<0.85))/length(priorchangeframe$uid)
percent[2]=100*sum((priorchangeframe$pA0.5>=0.85)*(priorchangeframe$pA0.5<=1))/length(priorchangeframe$uid)
thresholds=c("[0.5,0.85)","[0.85,1]")

thresholdsframe=data.frame(Theshold=thresholds,Percent=percent)

#check at each posterior how many people drop
#set the insignificant information threshold
#*************************
k=0
#**********************************
#set prior for checking

pri=0.85
infogatherframe0.85=data.frame(uid=vector('numeric'),shoulddrop=vector('numeric'),notdrop=vector('numeric'),kthreshdropped=vector('numeric'),betterthanA=vector('numeric'),significantchange=vector('numeric'))
for (i in 1:length(uidlist)){
pcfrow=subset(priorchangeframe, uid==uidlist[i])
playframe=subset(expframe,uid==uidlist[i]&priorA==pri)
sd=(pcfrow$pA0.5upper<pri)
#non significante version 
#sd=(pcfrow$pA0.5<pri)
nd=(pcfrow$pA0.5lower>pri)
#non significante version 
#nd=(pcfrow$pA0.5>pri)
ktd=(length(playframe$Aresponse)-sum(playframe$Aresponse)<=k)
btA=(sum(playframe$correct)>pri)
sigc=pcfrow$Prob0.5_0.85<0.05
addframe=data.frame(uid=uidlist[i],shoulddrop=sd,notdrop=nd,kthreshdropped=ktd,betterthanA=btA,significantchange=sigc)
infogatherframe0.85=rbind(infogatherframe0.85,addframe)
}

mu=c(0.85)
below0.85perc=100*sum(infogatherframe0.85$kthreshdropped*(infogatherframe0.85$shoulddrop))/sum(infogatherframe0.85$shoulddrop)
above0.85perc=100*sum(infogatherframe0.85$kthreshdropped*(infogatherframe0.85$notdrop))/sum(infogatherframe0.85$notdrop)

thresholdbelowmu=c(below0.85perc)
thresholdabovemu=c(above0.85perc)

droppingframe=data.frame(Mu=mu,Thresholdbelowmu=thresholdbelowmu,Thresholdabovemu=thresholdabovemu)
droppingframe=t(droppingframe)

#*********************************************
stargazer(droppingframe,colnames=FALSE,summary=FALSE)
#*********************************************

